In [15]:
%matplotlib inline
import numpy as np
import pandas as pd
import datetime
import matplotlib.pyplot as plt
import pylab
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
import statsmodels.api as sm
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf, q_stat
from statsmodels.graphics.gofplots import qqplot
from sklearn.metrics import mean_squared_error
from math import sqrt
In [16]:
np.random.seed(42)
data = pd.read_csv('../data/siteB.csv')
In [17]:
data['Date_Name'] = data['Date_Name'].apply(lambda x: datetime.datetime.strptime(x,'%Y-%m-%d'))
In [18]:
print data.info()
In [19]:
# Make a dataframe for each category, and put it into a list...len(list) = 16
# Essentially becomes a list of dataframe
def make_dataframe(df, lis):
f = lambda x: df[df['subCategory'] == x]
return map(f, lis)
temp_list = data['subCategory'].unique()
splitted_data = make_dataframe(data, temp_list)
print len(splitted_data)
print splitted_data[2][0:5]
In [20]:
# Looks like some have missing dates, like category D has only half of the dates.
idx = pd.date_range('09-26-2015', '09-30-2016')
for i, item in enumerate(temp_list):
splitted_data[i] = splitted_data[i].groupby(by=['Date_Name'])['sales'].sum()
print item, ':', splitted_data[i].shape
# We fill in the missing dates with value 0. Because this may be holiday and the shop made 0 sales.
splitted_data[i].index = pd.DatetimeIndex(splitted_data[i].index)
splitted_data[i] = splitted_data[i].reindex(idx, fill_value=0)
print item, ':', splitted_data[i].shape, '-', splitted_data[i].index
print '\n', splitted_data[0].head()
Step 2) Data exploratory¶
- We observe the shop has sales every single day, except for category D, which has sales for only 1/2 of the time.¶
- Weekly seasonal effect: weekend usually lower sales volume than weekday¶
- There is no upward or downward trend except for category E¶
- A convenient shop? The data is quite consistent throughout.¶
In [66]:
# Just to confirm my weekly theory is right,
# 2015-09-26, 2015-09-27, 2015-10-03, 2015-10-04, 2015-10-10, 2015-10-11, etc.
# are all weekends and have lower sales than weekdays.
# We will make a decomposition plot to further validate this.
print splitted_data[0].head(30)
In [65]:
# Time series in connected point plot format for each category.
# We can see the missing dates(or date with 0 sales) in the time series this way
fig, axs = plt.subplots(len(temp_list),1, figsize=(15, 96))
axs = axs.ravel()
for i, item in enumerate(data['subCategory'].unique()):
axs[i].plot(splitted_data[i].index, splitted_data[i], '-o')
axs[i].set_title('Category {} Sales'.format(item))
In [22]:
# Seasonal decomposition plots. Let's confirm the weekly seasonality.
def decomposition(timeseries, item):
decomposition = seasonal_decompose(timeseries, freq=7)
fig = plt.figure()
fig = decomposition.plot()
fig.set_size_inches(15, 8)
fig.suptitle('Category {} - Seasonal Decomposition Plot'.format(item))
_ = map(decomposition, splitted_data, temp_list)
Step 3.1) Stationarity Test¶
In [23]:
# We will use Dickey-Fuller Test to check the stationarity of our time series data.
def test_stationarity(timeseries, item):
#Determing rolling statistics
rol_mean = pd.rolling_mean(timeseries, window=12)
rol_std = pd.rolling_std(timeseries, window=12)
#Plot rolling statistics:
original = plt.plot(timeseries, color='blue',label='Original')
mean_ = plt.plot(rol_mean, color='red', label='Rolling Mean')
std_ = plt.plot(rol_std, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Category {} - Rolling Mean & Standard Deviation'.format(item))
plt.show(block=False)
#Perform Dickey-Fuller test:
print 'Results of Dickey-Fuller Test:'
adf_test = adfuller(timeseries, autolag='AIC')
adf_output = pd.Series(adf_test[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in adf_test[4].items():
adf_output['Critical Value (%s)'%key] = value
print adf_output
In [24]:
_ = map(test_stationarity, splitted_data, temp_list)
In [25]:
# We will use this function repeatedly to make plots for -> Transformed Data vs. Original Data
def plot_data(original, trans_2, item, title, label_1, label_2):
original.plot(label='{}'.format(label_1))
trans_2.plot(label='{}'.format(label_2))
plt.legend(loc='best')
plt.title('Category {} - {}'.format(item, title))
plt.show(block=False)
In [26]:
def log_data(timeseries, item):
# logged_data = timeseries.apply(lambda x: np.log(x))
return np.log(timeseries)
splitted_data_log = map(log_data, splitted_data, temp_list)
print splitted_data[0].head(), '\n', splitted_data_log[0].head()
In [27]:
# Let's see the log plot.
_ = map(plot_data, splitted_data, splitted_data_log, temp_list, \
['Log']*len(temp_list), ['original']*len(temp_list), ['logged']*len(temp_list))
In [28]:
# Replace inf, -inf with 0...because after we log(0), we got inf values
print splitted_data_log[14]['2015-12-20']
replace_inf = lambda x: x.replace([np.inf, -np.inf], 0)
splitted_data_log = map(replace_inf, splitted_data_log)
print splitted_data_log[14]['2015-12-20']
In [29]:
# Plot the log again after we replaced inf, -inf with 0.
_ = map(plot_data, splitted_data, splitted_data_log, temp_list, \
['Log']*len(temp_list), ['original']*len(temp_list), ['logged']*len(temp_list))
In [30]:
# Differencing to deal with the trend.
# Log first difference
def diff_data(timeseries, item):
timeseries_diff = timeseries - timeseries.shift()
return timeseries_diff
splitted_data_log_diff = map(diff_data, splitted_data_log, temp_list)
In [31]:
print splitted_data_log_diff[0].head()
In [32]:
_ = map(lambda x: x.dropna(inplace=True), splitted_data_log_diff)
_ = map(test_stationarity, splitted_data_log_diff, temp_list)
In [33]:
# Differencing to deal with seasonality...
# Log seasonal difference
splitted_data_log_diff_season = map(lambda x: x - x.shift(7), splitted_data_log_diff)
_ = map(lambda x: x.dropna(inplace=True), splitted_data_log_diff_season)
In [34]:
# Let's see the Dickey-Fuller test after taking the log, diff, & season diff...
_ = map(test_stationarity, splitted_data_log_diff_season, temp_list)
In [35]:
def acf_pacf_plot(timeseries, title, item, lags=30):
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(timeseries, lags=lags, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(timeseries, lags=lags, ax=ax2)
fig.suptitle('Category {} - {} ACF/PACF Plot'.format(item, title), fontsize=20)
In [37]:
# Observed data ACF/PACF plot.
_ = map(acf_pacf_plot, splitted_data, ['Observed Data']*len(temp_list), temp_list)
In [39]:
# Let's see the log transformed ACF/PACF plot.
_ = map(acf_pacf_plot, splitted_data_log, ['Log']*len(temp_list), temp_list)
In [41]:
# Diff log ACF/PACF plot.
_ = map(acf_pacf_plot, splitted_data_log_diff, ['Diff Log']*len(temp_list), temp_list)
In [43]:
# Finally, let's see the Season Diff Log transformed ACF/PACF plot.
_ = map(acf_pacf_plot, splitted_data_log_diff_season, ['Season Diff Log']*len(temp_list), temp_list, [range(0, 182, 7)]*len(temp_list))
In [47]:
# We will use this function to help us search the best p & q values, and seasonal P & Q values.
# Four for loops is a lot, but we will be using it mostly for 1 or 2 values in each p, q, P & Q values.
# Because we already have a good idea these values won't be more than 2.
def evaluate_models(dataset, p_values, q_values, P_values, Q_values):
best_score, best_cfg, best_season = float("inf"), None, None
for p in p_values:
for q in q_values:
for P in P_values:
for Q in Q_values:
order = (p,1,q)
season_order = (P, 1, Q, 7)
try:
mod = sm.tsa.statespace.SARIMAX(dataset, trend='n', order=order, \
seasonal_order=season_order)
result = mod.fit()
aic = result.aic
if aic < best_score:
best_score, best_cfg, best_season = aic, order, season_order
print('SARIMA%s %s AIC=%.3f' % (order, season_order, aic))
except:
continue
print('Best SARIMA%s %s AIC=%.3f' % (best_cfg, best_season, best_score))
return [best_cfg, best_season]
In [48]:
def model_building(timeseries, item, p=range(0,2), q=range(1, 2), P=range(0,2), Q=range(1, 2)):
orders = evaluate_models(timeseries, p, q, P, Q)
mod = sm.tsa.statespace.SARIMAX(timeseries, trend='n', order=orders[0], seasonal_order=orders[1])
results_1 = mod.fit()
print 'Category {}: '.format(item), results_1.summary()
return results_1
In [49]:
results_1 = map(model_building, splitted_data, temp_list)
In [50]:
# Next step: check the residual plots to see the models are good.
def get_residuals(timeseries_results):
residuals = pd.DataFrame(timeseries_results.resid)
return residuals
In [51]:
residuals = map(get_residuals, results_1)
In [52]:
def residual_analysis(residuals, item):
fig = plt.figure(figsize=(12,12))
ax1 = fig.add_subplot(411)
ax1 = sm.graphics.tsa.plot_acf(residuals, lags=30, ax=ax1, title='Category {} - ACF of Residual'.format(item))
ax2 = fig.add_subplot(412)
ax2.title.set_text('Category {} - Q-Q Plot of Std Residuals'.format(item))
ax2 = qqplot(residuals, line='r', ax=plt.gca())
# Plot residual errors.
ax3 = fig.add_subplot(413)
ax3.title.set_text('Category {} - Standardized Residuals'.format(item))
ax3 = plt.plot(residuals)
# residuals.plot(kind='kde')
temp_q_stat = sm.stats.diagnostic.acorr_ljungbox(residuals, lags=30)
ax4 = fig.add_subplot(414)
ax4.title.set_text('Category {} - P value for Ljung-Box Statistics'.format(item))
ax4 = plt.plot(temp_q_stat[1], 'o')
# print temp_q_stat[0]
plt.tight_layout()
In [53]:
_ = map(residual_analysis, residuals, temp_list)
The residuals plots mostly look good. Let's move on to forecasting.¶
In [ ]:
Step 5) Prediction In Sample, and 90 Days Ahead¶
- We will plot the forecast for each category¶
- Plot the relative contribution from each category¶
- Add them up and plot the total forecast¶
- Measure the RMSE¶
- We will also build a separate model and forecast using top-down approach for comparison (what we did so far is bottom up)¶
In [54]:
def plot_forecast(results, timeseries, item):
pred = (results.predict())
forecast_90 = (results.forecast(steps=90))
pred.plot(label='predict')
forecast_90.plot(label='forecast')
timeseries.plot(label='original')
plt.title(' Category {} Forecast '.format(item))
plt.legend()
plt.show(block=False)
return pred
In [55]:
pred = map(plot_forecast, results_1, splitted_data, temp_list)
In [56]:
# Report performance
def rmse_measurement(timeseries, pred, item):
mse = mean_squared_error(timeseries, pred)
rmse = sqrt(mse)
print('Category {} - RMSE: %.3f'.format(item) % rmse)
return rmse
In [57]:
rmse_total = map(rmse_measurement, splitted_data, pred, temp_list)
In [58]:
print 'Total RMSE: ', np.sum(rmse_total)
Total Forecast¶
In [59]:
forecast_90 = map(lambda x: x.forecast(steps=90), results_1)
total_forecast = pd.Series()
original_sales = pd.Series()
for i, df in enumerate(forecast_90):
if i == 0:
total_forecast = df
original_sales = splitted_data[i]
else:
total_forecast = total_forecast + df
original_sales = original_sales + splitted_data[i]
plt.plot(original_sales.index, original_sales, '-o')
plt.plot(total_forecast.index, total_forecast, '-o')
Out[59]:
In [60]:
def relative_cont_plot(forecast_90, item):
category = forecast_90 / total_forecast
# print 'Category {}: '.format(item), category.mean()
if category.mean() >= 0.08:
print 'Category {} Sales Contribution: '.format(item), category.mean()
_ = plt.plot(category.index, category, '-o', label=item)
_ = plt.legend(loc='right')
_ = plt.title('Relative Contribution from Each Category')
print 'Top contributors are...'
_ = map(relative_cont_plot, forecast_90, temp_list)
Top Down Approach - Build Model and Forecast on All Categories.¶
In [61]:
print original_sales.head()
idx = pd.date_range('09-26-2015', '09-30-2016')
original_sales.index = pd.DatetimeIndex(original_sales.index)
original_sales = original_sales.reindex(idx, fill_value=0)
print original_sales.shape, '-', original_sales.index
print '\n', original_sales.head()
In [ ]:
In [ ]:
In [79]:
total_forecast_1 = model_building(original_sales, 'Total', range(0, 3), range(0, 3), range(0, 3), range(0, 3))
In [80]:
total_pred_log = total_forecast_1.predict()
forecast_90_log = total_forecast_1.forecast(steps=90)
# total_pred_1.plot(label='pred')
original_sales.plot(label='original')
forecast_90_1.plot(label='forecast')
plt.legend()
rmse_measurement(original_sales, total_pred_1, 'Total')
Out[80]:
In [62]:
total_forecast = model_building(original_sales, 'Total')
In [63]:
total_pred = total_forecast.predict()
forecast_90 = total_forecast.forecast(steps=90)
total_pred.plot(label='pred')
original_sales.plot(label='original')
forecast_90.plot(label='forecast')
plt.legend()
Out[63]:
In [64]:
rmse_measurement(original_sales, total_pred, 'Total')
Out[64]:
In [ ]:
Summary¶
1) Statistical tests suggest¶
- We have lags and seasonal effects in the data¶
- So we need to use seasonal ARIMA model as opposed to regular ARIMA model¶
- The SARIMA models for each category will be slightly different, but shouldn't deviate too much from order 1¶
2) SARIMA models were used for forecasting each category¶
- Category A, D, Q: (0, 1, 1)x(1, 1, 1, 7)¶
- Category B, E, F, G, H, K, M, N, O: (1, 1, 1)x(0, 1, 1, 7)¶
- Category C, I, L: (1, 1, 1)x(1, 1, 1, 7)¶
3) We also used a top down approach (directly forecast without breaking down the categories)¶
- Top down approach (RMSE=2071)¶
-> Heavy grid search gave this model: (0, 1, 2)x(0, 1, 1, 7) RMSE=2073¶
-> Comparing to a simpler model (1, 1, 1, )x(0, 1, 1, 7) RMSE=2071, our theory of staying within the order 1 was correct!¶
- Bottom up approach (RMSE=2231)¶
- Seems like top down was able to generalize the time series better, it's interesting because one would expect the opposite¶
- But it could be the top down method is overfitting, so further validation is needed.¶
4) What could be improved?¶
- More data: if we have 5 years of data, we can definitely build more accurate models.¶
- More features: we could use weather data to see how climate affects sales in diff categories¶
- More features: also features like promotion or not, number of competitors nearby would also help¶
- Diff models: we could try xgboost¶
In [ ]: